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PREFACE 

The  work  reported  herein  was  conducted  at  Arnold  Engineering  Development  Center 
(AEDC),  Air  Force  Systems  Command  (AFSC).  The  work  was  accomplished  by  Sverdrup 
Technology,  Inc.,  AEDC  Group,  operating  contractor  for  the  propulsion  test  facilities  at 
AEDC,  AFSC,  Arnold  Air  Force  Base,  Tennessee,  under  Project  DC12EW.  The  Air  Force 
project  manager  was  K.  Zysk.  The  project  manager  for  Sverdrup  Technology  was  H.  T. 
Bentley,  III.  The  mansucript  was  submitted  for  publication  on  May  22,  1992. 
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1.0  INTRODUCTION 

Calculation  of  the  equilibrium  chemical  composition  of  complex  gas  mixtures  has  received 
much  attention  over  the  years.  Several  textbooks  (Refs.  1-6)  and  technical  reports/papers 
(Refs.  7-10)  present  the  necessary  fundamentals.  There  are  basically  two  approaches  toward 
computation  of  the  equilibrium  composition  for  a  complex  gas  mixture: 

1.  Equilibrium  constant  formulation, 

2.  Free-energy  minimization  formulation. 

Most  modern-day  multipurpose  chemical  equilibrium  computer  programs  (Refs.  1 1-14)  use 
the  free-energy  minimization  procedure.  These  programs  generally  employ  a  Newton-Raphson 
iteration  procedure  to  solve  a  set  of  simultaneous  nonlinear  algebraic  equations  with  various 
“practical”  enhancements  (such  as  restrictions  on  sizes  of  allowable  corrections  during  an 
iteration  cycle)  to  avoid  numerical  difficulties.  The  resulting  solution  minimizes  the  Gibbs 
free  energy  of  the  mixture  subject  to  the  constraint  of  conservation  of  elements.  A  linear 
equation  solver  /matrix  inversion  routine  must  be  invoked  at  each  iteration  cycle;  this  results 
in  most  of  the  computer  time  required  for  a  converged  solution  being  spent  to  perform  the 
requisite  numeric  matrix  operations. 

For  solution  of  complex  chemical  equilibrium  compositions  as  a  subroutine  in  a 
computational  fluid  dynamics  code  (say  an  unsteady  Navier-Stokes  flow-field  solver),  a  stable, 
fast  numerical  algorithm  which  does  not  require  matrix-type  numerical  operations  is  needed. 
This  becomes  especially  critical  when  it  is  realized  that  the  equilibrium  composition  must 
be  calculated  at  every  grid  point  (which  may  number  in  the  thousands)  over  numerous  time 
cycles  (which  may  also  number  in  the  thousands).  For  this  reason,  many  past  applications 
have  used  approximate  approaches  for  computing  complex  mixture  equilibrium  composition, 
the  work  by  Mascitti  (Ref.  15)  being  an  excellent  example  of  a  simplified  equilibrium 
hydrocarbon-air  combustion  model  specifically  created  for  use  in  air-breathing  engine  cycle 
computer  programs.  Other  examples  include  the  chemically  reacting  boundary-layer  programs 
by  Blottner  (Refs.  16-17)  which  contain  the  capability  for  chemical  equilibrium  ablation 
products  at  an  ablating  surface. 

The  various  programs  by  Spalding  and  associates  (Refs.  18-21)  for  analysis  of  hydrogen- 
air  combustion  under  various  flow  conditions  embody  a  significant  advance  in  the  numerical 
calculation  of  complex  equilibrium  composition  mixtures.  They  utilize  a  little-known  form 
of  the  quadratic  root  extraction  formula  to  develop  an  iterative  algorithm  for  solving  the 
hydrogen-air  equilibrium  composition  problem  based  upon  the  equilibrium  constant 
formulation  without  recourse  to  matrix  techniques.  This  approach  meets  the  desired  criteria 
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delineated  above  for  application  to  computational  fluid  dynamics  codes,  and  provides  the 
basis  for  the  current  work  which  formalizes  an  algorithmic  technique  that  is  computationally 
fast,  absolutely  stable,  and  simple  to  program  for  general  gas  mixtures.  Details  of  the  requisite 
quadratic  root  extraction  formula  are  presented  in  Appendix  A. 

2.0  APPROACH 

Calculation  of  the  equilibrium  composition  of  complex  gas  mixtures  using  the  equilibrium 
constant  formulation  generally  entails  the  following  five  steps  (see  Chap.  V,  Sec.  5  in  the 
book  by  Vincenti  and  Kruger,  Ref.  4): 

1.  Specify  an  independent  set  of  chemical  reactions  containing  all  the  chemical 
species  of  interest  for  the  gas  mixture. 

2.  Apply  the  law  of  mass  action  individually  to  each  of  the  chemical  reactions 
from  Step  1. 

3 .  Supplement  the  equations  from  Step  2  by  equations  expressing  the  conservation 
of  atomic  nuclei  and  the  fact  of  zero  net  electronic  charge  if  ionization  is  present. 

4.  Combine  the  equations  resulting  from  Steps  2  and  3  to  form  a  set  of 
simultaneous  nonlinear  algebraic  equations,  the  number  of  equations  being  equal 
to  the  number  of  chemical  species  in  the  gas  mixture. 

5.  For  a  specified  pressure,  temperature,  and  initial  composition,  solve  the 
equations  from  Step  4  to  yield  the  gas  mixture  equilibrium  composition. 

It  is  only  in  Step  5  that  the  actual  numerical  solution  algorithm  becomes  important.  However, 
the  manner  in  which  the  governing  equations  are  posed  in  Step  4  is  critical  to  the  selection 
and  application  of  this  solution  algorithm. 

As  discussed  in  Section  1.0,  the  present  work  is  concerned  with  a  fast,  absolutely  stable 
numerical  algorithm  for  solving  the  equilibrium  composition  problem  without  recourse  to 
matrix  techniques.  This  approach  will  be  illustrated  through  application  to  two  complex  gas 
mixtures: 

1.  High-temperature  ionized  air, 

2.  Hydrocarbon-air  products  of  combustion. 

Sufficient  detail  will  be  included  to  allow  application  of  the  algorithm  to  other  mixtures. 
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3.0  APPLICATION  EXAMPLES 

3.1  HIGH-TEMPERATURE  IONIZED  AIR 

Following  Ref.  4,  ionized  air  at  temperatures  up  to  about  8,000  K  may  be  regarded  as 
composed  of  seven  chemical  species:  O2,  O,  N2,  N,  NO,  NO+ ,  and  e  — .  The  equilibrium 
of  the  resulting  mixture  can  be  taken  to  be  controlled  by  the  following  four  independent 
chemical  reactions,  identified  by  the  letter  r: 

r  =  1:  02**20 
2:  N2^2N 
3:  NO**N  +  O 
4:  NCMNO+  +  e~ 

Application  of  the  law  of  mass  action  to  these  four  reactions  yields  the  conditions  for 
equilibrium  of  the  mixture  as 


r 


Xo  =I„ 
Xq2  G  Kc' 


(1) 


2: 


Xti  1 


C) 


3: 


XNX0  _  1 
XnO  e  Cj 


(3) 


4: 


Xnq+  Xg- 

XNO 


(4) 


where  Xi  denotes  the  number  of  moles  of  the  species  i  per  unit  mass  of  the  mixture,  q  denotes 
the  density  of  the  mixture,  and  K^.r  is  the  equilibrium  constant  for  the  reaction  r  based  on 
concentration.  These  equations,  which  are  four  equations  in  seven  unknowns,  must  be 
supplemented  by  three  additional  equations  expressing  the  atomic  conservation  of  oxygen 
and  nitrogen  nuclei  as  well  as  the  fact  of  zero  net  electronic  charge.  These  are,  in  terms  of 
the  X  mote-mass  ratios: 


2Xq2  +  Xq  +  XNc>  +  XNO+  -  2Xq2 


(5) 
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2XNz  +  XN  +  XNO  +  XNQ+  =  2Xn2  (6) 

^N0+  —  ^e“  =  ®  (?) 

where  Xo^and  XN’  are  the  02  and  N2  compositions  at  low  temperature  where  the  mixture 
is  totally  02  and  N2.  Values  for  Xq2  and  XN2  that  are  a  good  approximation  to  real  air  are 


Xo2  =  0.00733 


mole 

gram  -  mixture 


0.0273 


mole 

gram  -  mixture 


as  given  in  Ref.  4.  With  the  pressure  p  and  temperature  T  of  the  mixture  specified,  the  mixture 
density  q  is  given  from  the  thermal  equation  of  state  for  the  gas  mixture 


_  P 
e  E  XjRT 

i 


(8) 


where  R  is  the  universal  gas  constant  and  the  summation  is  over  the  gaseous  species  of  the 
mixture. 


Substitution  of  Eqs.  (1),  (2),  and  (3)  allows  Eqs.  (5)  and  (6)  to  be  written  as 


r  2q  i 

•  -2 

eXN- 

.  1 

[kJ 

x0  + 

1  +  K 
kc3J 

x0  + 

2Xo2  +  ^no-J 

r  2e  1 

-.2 

r  eXoi 

Xn  + 

.  1 

[kJ 

Xn  + 

1  +  K 
KC3J 

2Xn2  +  XNO+j 

(9) 


(10) 


which  are  of  the  proper  form  for  application  of  the  quadratic  root  extraction  discussed  in 
Appendix  A.  Note  that  Eqs.  (9)  and  (10)  result  in  A,  B,  and  C  coefficients  which  satisfy 
the  necessary  criteria  discussed  in  the  Appendix;  i.e.,  A  and  B  are  both  positive  and  C  is 
negative  (it  is  physically  impossible  for  Xno+  to  be  larger  in  magnitude  than  Xq2  or  Xn2). 
A  pseudo-coding  of  a  suitable  algorithm  for  numerical  solution  of  these  equations  is  presented 
in  List  I .  Note  that  the  algorithm  is  very  simple  (successive  substitution  with  iteration  to 
convergence).  The  computational  philosophy  embodied  in  the  algorithm  is  detailed  in  four 
steps: 


1.  Solve  first  for  the  atomic  species  compositions  from  the  quadratic  equations 
which  result  through  substitution  of  the  law  of  mass  action  equations  into  the 
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conservation  of  atomic  nuclei  equations.  Use  the  quadratic  root  extraction 
formula  of  Appendix  A. 

2.  Next,  solve  for  the  diatomic  species  compositions  by  substitution  of  the  atomic 
species  compositions  from  Step  1  into  the  law  of  mass  action  equations. 

3 .  Finally,  update  the  density  of  the  gas  mixture  based  upon  the  thermal  equation 
of  state. 

4.  Iterate  to  convergence  on  all  species  compositions. 

Also  illustrated  in  the  algorithm  pseudo-coding  of  List  1  is  the  necessary  specification  of 
the  initial  composition,  as  well  as  the  curve  fit  form  of  the  equilibrium  constants  taken  from 
Ref.  4  for  the  four  chemical  reactions  in  question. 

Practical  application  of  this  algorithm  reveals  convergence  within  three  to  five  iteration 
cycles  using  reasonable  initial  guesses  for  the  species  mole-mass  ratios  and  gas  density.  For 
poor  initial  guesses  (zero  values  for  the  mole-mass  ratios  and  gas  density  based  on  the  low- 
temperature  O2  -  N2  air  mixture),  convergence  is  obtained  in  from  four  to  eight  iteration 
cycles. 

3.2  HYDROCARBON-AIR  PRODUCTS  OF  COMBUSTION 

Many  combustion  problems  require  the  knowledge  of  the  equilibrium  composition  of 
a  complex  hydrocarbon-air  mixture  which  may  contain  water  vapor  and  various  oxides.  For 
the  present  example  the  hydrocarbon-air  mixture  will  be  taken  to  be  composed  of  the  following 
thirteen  chemical  species:  C,  H,  N,  O,  O2,  N2,  H2,  OH,  H20,  CO,  C02,  NO,  and  N02.  The 
equilibrium  of  the  resulting  mixture  can  be  taken  to  be  controlled  by  the  following  nine 
independent  chemical  reactions,  identified  by  the  letter  r: 


r  =  1:  Oa  =  O  +  O 


2:  N2  =  N  +  N 


3:  H2  =  H  +  H 


4:  OH  =  O  +  H 


5:  H20  =  H  +  H  +  O 


9 


AEDC-TR-91  -14 


6:  CO  =  C  +  O 
7:  C02  =  C  +  O  +  O 
8:  NO  =  N  +  O 


9:  N02  =  N  +  O  +  O 

Applying  the  law  of  mass  action  to  these  eight  reactions  yields  the  conditions  for  equilibrium 
of  the  mixture  as 


,  Xo  1 
r  =  1:  77-  -  -  KCl 
X02  Q 


01) 


Xn  1 

2:  ^  =  "Kc2 
Xn2  e  i 


(12) 


Xh  1 

3:  xiL"IKc* 
xh2  c 


(13) 


4:  **5.1*. 

Xqh  G  4 


5.  ^2  =  1 

XH2o  G2  Kcs 


(14) 

(15) 


6:  ^  =  ^KC6 
Xco  G  6 


(16) 


xcx6  1 

7:  =  -3  Kc, 

XC02  g2  Ct 


(17) 


8:  **2.IK 

Xno 


C8 


(18) 


Xnx6_1 

XNq2  C2  KC9 


(19) 
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where  X[  denotes  the  number  of  moles  of  the  species  i  per  unit  mass  of  the  mixture,  e  denotes 
the  density  of  the  mixture,  and  K<;r  is  the  equilibrium  constant  for  the  reaction  r  based  on 
concentration.  These  equations,  which  are  nine  equations  in  thirteen  unknowns,  must  be 
supplemented  by  four  additional  equations  expressing  the  atomic  conservation  of  oxygen, 
nitrogen,  hydrogen,  and  carbon  nuclei.  These  are,  in  terms  of  the  X  mole-mass  ratios: 


2Xo2  +  Xo  4  XH2o  +  Xoh  +  Xno  +  2Xno2  4  2Xcq2  +  Xco  =  Xq  (20) 

2XN2  +  Xn  +  Xn0  +  XN02  =  Xn  (21) 

2XH2  +  Xh  4  2XH2o  4  Xoh  -  X*„  (22) 

Xc  4  Xco2  4  Xco  =  Xc  (23) 


where  Xq,  X^,  X^,  and  X^  represent  the  initial  elemental  composition  of  the  mixture. 

At  this  point  exactly  the  same  computational  philosophy  discussed  earlier  relative  to  the 
ionized  air  example  is  followed,  namely  substitution  of  the  law  of  mass  action  equations 
into  the  four  conservation  of  atomic  nuclei  equations  to  yield  three  quadratic  equations  for 
the  atomic  oxygen  O,  atomic  nitrogen  N,  and  atomic  hydrogen  H  mole-mass  ratios,  and 
a  single  equation  for  the  atomic  carbon  C  mole-mass  ratio.  These  three  quadratic  equations 
are  solved  using  the  quadratic  root  extraction  formula  of  Appendix  A  with  the  results 
substituted  into  the  law  of  mass  action  equations  to  yield  the  diatomic  H2,  02,  N2,  OH,  NO, 
and  triatomic  H20,  C02,  N02  mole-mass  ratios.  A  pseudo-coding  of  this  algorithm  is  given 
in  List  2;  note  the  similarity  with  the  ionized  air  algorithm  of  List  1.  Iterative  convergence 
of  this  algorithm  to  the  final  mixture  composition  can  be  significantly  enhanced  by: 

1.  Performing  an  inner  iteration  within  the  overall  iteration  process  for  some 
specified  number  of  cycles  (typically  five)  on  the  four  conservation  of  atomic 
nuclei  equations. 

2.  Extrapolating  the  atomic  oxygen  O  mole-mass  ratio  based  upon  the  two  most 
recent  iterations. 

These  two  “practical”  considerations  are  embodied  in  the  algorithm  formulation. 
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Application  of  the  present  successive  substitution  algorithm  to  numerous  hydrocarbon- 
air  product  of  combustion  mixtures  for  various  fuels  (e.g.,  methane,  hydrogen,  kerosene, 
and  JP-4)  reveals  convergence  within  three  to  five  iterations  using  reasonable  initial  guesses 
for  the  species  mole-mass  ratios  and  gas  density.  For  poor  initial  guesses  (zero  values  for 
the  mole-mass  ratios  and  gas  density  based  on  the  low-temperature  O2  -N2  air  mixture  only), 
convergence  is  obtained  in  from  four  to  eight  iteration  cycles.  The  algorithm  has  been  exercised 
over  a  wide  range  of  fuel-air  ratios  (fuel  lean  to  stoichiometric  to  fuel  rich)  with  no  degradation 
in  performance.  In  general,  it  has  proven  to  be  highly  reliable  and  robust,  mainly  because 
of  the  desirable  numerical  properties  of  the  quadratic  root  extraction  formula  of  Appendix  A. 

4.0  ALGORITHM  VALIDATION 

To  validate  the  presently  proposed  algorithm  for  application  to  complex  gas  mixtures 
equilibrium  compositions,  calculations  have  been  performed  for  two  hydrocarbon-air  mixtures 
using  (1)  the  hydrocarbon-air  iterative  successive  substitution  algorithm  of  List  2,  (2)  the 
free-energy  minimization  procedure  of  NASA  SP-273  (Ref.  13)  involving  matrix  inversion, 
and  (3)  the  equilibrium  constant  approach  of  AEDC-TR-7 1-256  (Ref.  14)  involving 
linearization  of  the  governing  conservation  equations  in  conjunction  with  Newton-Raphson 
iteration.  All  three  approaches  utilized  exactly  the  same  set  of  thermochemical  properties, 
namely  the  JANNAF  (Joint  Army,  Navy,  NASA,  Air  Force)  data  curve  fits  from  Ref.  13. 
Hence,  any  differences  observed  upon  comparison  of  results  from  the  three  totally  different 
methods  can  be  directly  attributed  to  numerics  and/or  precision  in  performing  the  actual 
calculations. 


Mixture 

Pressure,  atm 

Temperature,  K 

Theoretical  Air, 
percent 

Methane/Air 

1.0 

1,000.0 

110.0 

Propane/Air 

50.0 

2,500.0 

100.0 

As  shown  in  Table  1  for  the  methane/air  mixture,  the  present  successive  substitution 
hydrocarbon-air  products  of  combustion  algorithm  yields  results  in  perfect  agreement  with 
the  other  two  methods.  Only  four  iterations  were  required  for  convergence  of  the  present 
algorithm  with  an  initial  guess  of  zero  for  all  species  mole-mass  ratios  and  the  mixture  density 
initialized  at  the  low-temperature  02-N2  air  value.  Table  2  presents  corresponding  results 
for  the  propane/air  mixture;  note  that  the  disagreement  between  any  of  the  three  methods 
for  any  species  is  in  the  fifth  decimal  digit.  For  this  propane/air  case,  five  iterations  were 
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required  for  convergence  of  the  present  algorithm  using  initial  guesses  as  described  above 
for  the  methane/air  mixture. 

As  an  additional  validation,  the  presently  proposed  algorithm  was  used  to  calculate  the 
adiabatic  flame  temperature  and  resulting  chemical  composition  for  a  gaseous  rocket 
combustion  chamber  application;  specifically: 

Fuel:  75-percent  ethylene  and  25-percent  nitrogen  (by  volume) 

Oxidizer:  Oxygen 

Oxidizer /Fuel  Ratio:  2.14 

Chamber  Pressure:  6.8  atm 

using  the  hydrocarbon-air  iterative  successive  substitution  algorithm  and  the  free-energy 
minimization  procedure  of  NASA  SP-273  involving  matrix  inversion.  Both  approaches  used 
exactly  the  same  thirteen  chemical  species  and  same  set  of  thermochemical  properties,  namely 
the  JANNAF  data  curve  fits  from  Ref.  13.  Hence,  any  differences  observed  upon  comparison 
of  results  from  the  two  totally  different  methods  can  be  directly  attributed  to  numerics  and/or 
precision  in  performing  the  actual  calculations. 

Table  3  shows  the  present  hydrocarbon-air  products  of  combustion  algorithm  yields  results 
in  perfect  agreement  with  the  NASA  SP-273  chamber  calculation  procedure.  For  a  given 
value  of  the  adiabatic  flame  temperature,  three  to  four  iterations  were  required  for  convergence 
of  the  present  algorithm  using  the  species  mole-mass  ratios  from  the  previous  adiabatic  flame 
temperature  calculation  as  the  initial  guess.  Six  iterations  on  the  adiabatic  flame  temperature 
were  required  for  adiabatic  flame  convergence  starting  from  an  initial  guess  of  3,000  K.  Nine 
iterations  in  the  NASA  SP-273  chamber  calculation  procedure  were  required  for  adiabatic 
flame  convergence. 

The  above  results  clearly  show  that  the  presently  proposed  interactive  successive  substitution 
algorithm  does,  indeed,  converge  rapidly  to  the  correct  solution  for  complex  hydrocarbon 
mixtures,  even  with  poor  initial  guesses  to  start  the  iteration  process.  Better  initial  guesses 
significantly  increase  the  rate  of  iterative  convergence  as  discussed  previously. 

Execution  timing  studies  of  the  presently  proposed  iterative  successive  substitution 
algorithm  given  in  List  2  relative  to  the  NASA  SP-273  algorithm  reveal  the  following 
comparison: 
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Relative  Time  for 

Algorithm  One  Iteration 


Present  Successive  1.0 

Substitution  Algorithm 

NASA  SP-273  13.2 

for  the  gaseous  rocket  combustion  chamber  application  discussed  earlier.  Both  algorithms 
used  the  exact  same  number  of  chemical  species  (namely  thirteen),  and  the  relative  time  for 
one  iteration  was  determined  via  system  time  calls  at  the  start  and  finish  of  each  iteration 
cycle.  Only  the  time  spent  in  actually  performing  the  equilibrium  composition  numerics  was 
included;  no  input/output  or  thermodynamic  property  calculations  were  contained  in  the 
timing.  Thus,  the  presently  proposed  algorithm  is  over  an  order  of  magnitude  faster  than 
the  NASA  SP-273  algorithm  (which  requires  matrix  inversion)  for  this  thirteen-species  case. 

For  application  of  the  present  successive  substitution  algorithm  in  computational  fluid 
dynamic  codes,  a  simple,  easy-to-use  curve-fit  format  for  the  equilibrium  constant  Kc  is 
desired.  Appendix  B  presents  an  Arrhenius-type  curve  fit  (Ref.  4,  pp.  214-215)  to  the  JANNAF 
data  curve  fit  from  Ref.  13.  As  shown  in  Appendix  B,  the  Arrhenius-type  curve  fit  is  an 
excellent  approximation  to  the  JANNAF  data  curve  fit  for  all  chemical  species/reactions 
considered.  When  they  are  applied  to  the  present  algorithm  validation  cases  of  Tables  1  and 
2,  the  Arrhenius-type  curve  fit  results  are  in  essentially  identical  agreement  with  the  JANNAF 
data  curve  fit  values  from  Tables  1  and  2  as  shown  in  Table  4.  Hence,  the  equilibrium  constant 
curve  fit  coefficients  given  in  Table  B-l  are  recommended  for  general  use  where  many 
evaluations  of  the  equilibrium  constant  Kc  must  be  performed,  as  in  a  computational  fluid 
dynamics  code  at  numerous  grid  points  over  numerous  time  steps. 

5.0  CONCLUDING  REMARKS 

A  simple  numerical  algorithm  has  been  developed  for  the  computation  of  the  equilibrium 
composition  of  complex  gas  mixtures  based  upon  application  of  the  equilibrium  constant 
approach.  The  algorithm  utilizes  a  little-known  form  of  the  quadratic  root  extraction  formula 
to  solve  for  the  gas  mixture  equilibrium  composition  using  a  successive  substitution  with 
iteration  to  convergence  technique.  Experience  with  the  algorithm  for  various  complex  gas 
mixtures  (ionized  air  and  hydrocarbon  products  of  combustion)  has  revealed  excellent  stability 
and  convergence  properties.  Its  application  for  arbitrary  gas  mixtures  in  computational  fluid 
dynamics  codes  as  an  economical,  easy-to-use  subroutine  is  straightforward. 
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List  1.  Algorithm  Pseudo-Coding  for  High-Temperature  Ionized  Air 


Equilibrium  Composition  of  High  Temperature  Ionized 
Chemical  Species ■  O,  02.  N.  M2,  NO,  NO* ,  e- 

Air 

Chemical  Reactions: 

r 

02  -  0  *  O 

i 

N2  -  N  «■  N 

2 

HO  -  N  +  O 

3 

HO  -  NO+  +■  a- 

4 

(Lav  Temperature  Mole-Mase  Ratio  Composition  of  D2  -  N2  Air) 

XSTAR  f  02  |  •  0.00723 
XSTAR  |  N2  |  =  0.0273 

(Calculate  Equilibrium  Constants  For  Each  Of  The  Four  Chemical  Reactions) 

KC  [  1  ]  -  1200.0  *  SXP<  -59500 , 0  /  T  )  /  SQRT(  T  ) 

KC  [  2  )  *  18. 0  »  EXP(  -113000.0  /  T  ) 

KC  (  3  )  =■  4.0  •  BXP(  -75500.0  /  T  ) 

KC  (  4  I  =»  14.4E-1Q  *  T  *  SQRT(  T  J  *  SXP(  -107000.0  /  T  ) 

(Repeat  Iterative  Solution  Procedure  For  Hole-Mass  Ratios  Until  Converged) 

8SGIH  REPEAT  UNTIL  CONVERGED  ON  ALL  X  (  I  J  FOR  ALL  SPECIES 

(Solve  Conservation  of  Atonic  Element  Equations) 

A  a  2.0  •  RHO  /  KC  C  1  ] 

B  =  1.0  +  RHO  »  X  (  V  ]  /  KC  (  3  ) 

C  >  -2.0  *  XSTAR  I  02  ]  +  X  I  NO+  ] 

X  [  0  ]  -  -2.0  «  C  /  (  B  +  SQRT(  B  *  B  -  4.0  *  A  *  C  )  ) 

A  a  2.0  *  RBO  /  KC  [  2  ] 

a  -  1.0  +  RHO  *  X  I  0  )  /  KC  (  3  ] 

C  »  -2.0  *  XSTAR  [  N2  ]  +  X  [  NO+  ] 

X  [  N  ]  «  -2.0  *  C  /  {  B  +  SORT |  8  •  B  -  4.0  •  A  *  C  )  ) 

(Solve  Law  Of  Mass  action  Algebraic  Equations} 

X  [  02  i  =  RHO  *  X  [  0  I  *  X  [  0  ]  /  KC  (  1  ] 

X  [  N2  |  -  RHO  •X(N|*X[N]/KC[2] 

X  [  NO  I  =  HHO  *X(N)*X(0)/KC[3) 

X  [  NO+  ]  *  SORT {  KC  (  4  ]  *  K  (  NO  1  /  RHO  ) 

{  Apply  Conservation  of  Electronic  Charge  } 

X  (  e-  )  »  X  [  NO*  ] 

(Calculate  Density  Of  Gas  Mixture  Using  The  Thermal  Equation  of  State) 

SUM  >0.0 

BEGIN  DO  LOOP  FOR  I  =  1  TO  NUMBER  OF  GASEOUS  SPECIES 
SUM  -  SUM  +  X  [  I  1 

END  DO  LOOP  FOR  I  >  1  TO  NUMBER  OF  GASEOUS  SPECIES 
RHO  -  P  /  (  SUM  *  R  •  T  J 

(insert  Suitable  Logic  For  Testing  Convergence  of  Iteration  Process  Here) 
BND  REPEAT  UNTIL  CONVERGED  ON  ALL  X  [  I  J  FOR  ALL  SPECIES 
(  End  of  Algorithm  ) 
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List  2.  Algorithm  Pseudo-Coding  for  Hydrocarbon- Air  Products  of  Combustion 


Equilibrium  Composition  For  C-H-N-O  Mixtures 
Hydrocarbon  -  Air  Products  of  combustion 
Chemical  Species i  C,  H,  H,  0,  02.  Hi.  K2.  OH.  BIO,  CO,  C02,  HO.  HO 2 

Chemical  React  Ions i  r 

Ol  »  O  ♦  0  1 

N2  a  N  *  N  3 

H2  -  H  +  H  3 

OH  -  O  +  H  4 

H20  -  H  +  H  *  O  S 

co  ■  c  t  o  a 

C02  »  C  ♦  O  ♦  O  7 

HO  -  M  ♦  0  a 

_ N02  a  H  ‘  O  ♦  O  9 


(Repeat  Iterative  Solution  Procedure  For  Mole-Mass  Ratios  Until  converged) 
BEGIN  REPEAT  UNTIL  CONVERGED  ON  ALL  X  J  I  |  FOR  ALL  SPECIES 
(Solve  Conservation  Of  Atomic  Element  Equations) 

BEGIN  DO  LOOP  PGR  J  ■  1  TO  HUMBER  OP  INNER  ITERATIONS 
A  “  2.0  *  ABO  *  |  1.0  /  KC  (  1  1 

MUffl  •  (  *  I  C  )  /  KC  [  7  ]  t  *  (  H  )  /  «c  (  9  I  I  ) 

H  -  1.0  +  RHO  *  <  X  [  H  ]  /  KC  (  4  ]  +  X  (  C  ]  /  KC  [  S  ] 

♦  I  [  »  ]  /  K  [  8  |  *  RKO  *X[B)*X(H]/KC(Sj) 

C  a  — XSTAR  (0| 

100  a  XO 
XO  -  X  [  0  ] 

X(O]--2.0*C/(8+  SQRTl  B‘B-4.0*A*C)  ( 

CAPFP-  CAFF 

CAPF  -  XO  -  X  [  0  J 

IP  (  (  J  >  2  1  AND  (  CAPF  <>  CAPPF  )  )  THEN 

EXTRAPOLATE  -  XO  -  CAPF  »  (  XO  -  XOO  J  /  (  CAPF  -  CAPPF  > 

IP  (  EXTRAPOLATE  >  0.0  )  THEN  I  |  0  I  •  EXTRAPOLATE 
END  If  THEN 


A  -  2.0  *  RHO/  KC  (  2  1 

B  *  1.0  +  RHO  *  X  [  O  )  »  (  1.0  /  KC  [  B  ]  +  RHO  •  X  (  0  )  /  KC  [  9  )  I 
C  ■  -XSTAR  (  N  ) 

X  [  N  1  ■  -2.0  *  C  /  (  B  *■  SORT!  B»B-4.0»A»C)  > 


A  -  2.0  *  RHO  *  (  1.0  /  KC  [  3  J 
B  -  1.0  ♦  RB0  *  X  (  0  J  /  KC  (  4 
C  -  -XSTAR  [  H  ] 

X  (  B  |  ■  -2,0  *  C  /  (  B  *  S0RT( 


+RHO*X[OJ/KC[S1) 
B»fl-4.0»A*C)  ] 


X  [  C  |  -  XSTAR  [  C  J  /  |  1.0  «■  RBO  *  X  (  O  J  »  (  1.0  /  BQKC  [  6  ] 
+  RHO*X[0)/KC[7]J) 


END  DO  LOOP  FOR  J  a  1  TO  NUMBER  OF  INNER  ITERATIONS 


(Solve  Law  of  Haas  Action  Algebraic  Equations) 


X  [  02  I  «  RHO  *X(0]  1  I  (  O  )  /  RC  | 
X  [  N2  |  a  R|£0  *  X  [  N  ]  a  X  [  N  |  /  XC  [ 
X  [  B2  ]  »  RHO  *X|K]»X(H]/KC( 
X(OK]aRHO*X|0]*X[Hl/KC| 
X  [  920  J  »  RHO  »  RBO  *  X  (  H  ]  •  X  I  K  ] 
*  [  CO  |  ■  RBO  *  I  |  C  ]  '  I  [  O  |  /  KC  ( 
X  (  C02  ]  -  RHO  *  RHO  *  X  (  C  ]  *  X  (  O  ) 
I  [  NO  |  =  RHO  *  X  (  N  ]  *  X  (  0  )  /  KC  [ 
X  (  N02  ]  -  RBO  *  RHO  *X[N]*XIO] 


(  O 
[  O 
(  0 


1  /  KC  [  5  1 
]  /  KC  t  7  ) 
I  /  KC  (  9  ] 


(Calculate  Density  Of  Gas  Mixture  Using  The  Thermal  Equation  of  State) 


SUM  -0.0 

BEGIN  DO  LOOP  FOR  I  ■  1  TO  NUMBER  OF  GASEOUS  SPECIES 
SUM  “  SUM  ♦  X  (  I  J 

END  DO  LOOP  FOR  I  a  I  TO  NUMBER  OF  GASEOUS  SPECIES 
RHO  "  F  /  (  SUM  «  R  »  T  ) 

(insert  Suitable  Logic  For  Testing  Convergence  of  Iteration  Process  Here) 
END  REPEAT  UNTIL  CONVERGED  OF  ALL  X  (  I  ]  FOR  ALL  SPECIES 
(End  Of  Algorithm) 
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Table  1.  Equilibrium  Composition  of  Methane/Air  Mixture* 


Species 

(Mass  Fraction) 

Successive 

Substitution 

Algorithm 

NASA  SP-273 
(Ref.  13) 

AEDC-TR-7 1-256 
(Ref.  14) 

02 

0.02010 

0.02010 

0.02010 

h2o 

0.11320 

0.11320 

0.11320 

n2 

0.72841 

0.72841 

0.72841 

NO 

0.00001 

0.00001 

0.00001 

co2 

0.13827  0.13827 

0.13827 

C,  0,  H2  H,  OH,  N,  N02,  CO 

<  0.00001 

<  0.00001 

<  0.00001 

*  Pressure  =  1  aim 
Temperature  =  1,000  K 
110-percent  Theoretical  Air 

Air  Taken  to  be  21-percent  O2  and  79-percent  N2  by  Volume 
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Table  2.  Equilibrium  Composition  of  Propane/Air  Mixture* 


Species 

(Mass  Fraction) 

Successive 

Substitution 

Algorithm 

NASA  SP-273 
(Ref.  13) 

AEDC-TR-7 1-256 
(Ref.  14) 

02 

0.00404 

0.00402 

O 

0.00007 

0.00007 

h2 

0.00015 

0.00015 

0.00015 

h2o 

0.09638 

0.09640 

0.09638 

OH 

0.00148 

0.00143 

0.00148 

N2 

0.71932 

0.71931 

0.71933 

NO 

0.00319 

0.00320 

0.00319 

C02 

0.16622 

0.16624 

0.16621 

CO 

0.00916 

0.00915 

0.00916 

C,  H,  N,  N02 

<  0.00001 

<  0.00001 

<  0.00001 

*  Pressure  =  50  atm 
Temperature  =  2,500  K 
100- percent  Theoretical  Air 

Air  Taken  to  be  21-percent  O2  and  79-percent  N2  by  Volume 
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Table  3.  Adiabatic  Flame  Temperature  Composition  for  a  Gaseous  Rocket  Combustor* 


Successive 

Substitution 

Algorithm 

NASA  SP-273 
(Ref.  13) 

Adiabatic  Flame  Temperature  (K) 

3431.9 

3431.9 

Species  (Mole  Fraction) 

H 

0.0518 

0.0518 

O 

0.0339 

0.0339 

o2 

0.0463 

0.0463 

n2 

0.0602 

0.0602 

h2 

0.0726 

0.0726 

OH 

0.0819 

0.0820 

h2o 

0.2519 

0.2519 

CO 

0.2727 

0.2726 

co2 

0.1186 

0.1187 

NO 

0.0101 

0.0101 

C,  N,  N02 

<  0.0001 

<  0.0001 

*  Pressure  =  6.8  atm 

Oxidizer/Fuel  Ratio  =  2.14 

Fuel:  75-percent  Ethylene  and  25-percent  Nitrogen  (by  Volume) 
Oxidizer:  Oxygen 
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Table  4.  Influence  of  Equilibrium  Constant  on  Equilibrium  Composition 
Using  Successive  Substitution  Algorithm 


Species 

(Mass  Fraction) 

Methane/Air  Mixture* 

Propane/Air  Mixture* 

JANNAF  Data  Kc 

Appendix  B  Kc 

JANNAF  Data  IQ 

Appendix  B  Kc 

o2 

0.02010 

0.02010 

0.00402 

0.00402 

o 

<  0.00001 

<  0.0000 1 

0.00007 

0.00007 

Hz 

<  0.0000 l 

<  0.00001 

0.00015 

0.00015 

h2o 

0.11320 

0.11320 

0.09638 

0.09639 

OH 

<  0.00001 

<  0.00001 

0.00148 

0.00147 

N2 

0.72841 

0.72841 

0.71932 

0.71932 

NO 

0.00001 

0.00001 

0.00319 

0.00319 

co2 

0.13827 

0.13827 

0.16622 

0.16623 

CO 

<  0.00001 

<  0-.00001 

0.00916 

0.00915 

C,  H,  N,  NOj 

<  0.00001 

<  0.00001 

<0.00001 

<  0.00001 

*  Pressure  -  1  atm 
Temperature  =  1,000  K 
110-percent  Theoreiical  Air 
Air  Taken  to  be  21-percent  02  and 
79-percent  N2  by  Volume 


+  Pressure  =  50  atm 
Temperature  =  2,500  K 
100-percent  Theoreiical  Air 
Air  Taken  to  be  21-percent  02  and 
79-percent  N2  by  Volume 
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APPENDIX  A 

QUADRATIC  ROOT  EXTRACTION  FORMULA 

The  two  real  roots  of  a  quadratic  equation 

AX2  +  BX  +  C  =  0 


are  well  known  to  be 

-  B  ±  -s/b2  -  4 AC 


What  is  not  so  well  known  is  that  these  two  real  roots  are  also  given  by* 

-2C 

X  =  - - 

B  ±  V®2  “  4AC 

which  has  many  desirable  advantages  for  numerical  applications.  If  it  is  physically  known 
that  both  A  and  B  are  greater  than  zero  (i.e.,  positive)  and  that  C  is  less  than  zero  (i.e., 
negative),  then  the  positive  root  of  the  equation  (which  is  the  root  of  interest  in  most  physical 
applications)  is  given  by 


-2C 

X  =  - . 

B  +  VB2  -  4AC 

without  any  danger  of  numerical  precision  problems  caused  by  subtraction  error.  Further 
observe  that  this  form  reduces  to  the  correct  positive  root  in  the  limit  of  either  the  A  or  B 
coefficient  approaching  zero. 


*A1AA  Aerospace  Design  Engineers  Guide,  Revised  and  Enlarged,  January  1987,  pp.  1-3. 
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APPENDIX  B 

EQUILIBRIUM  CONSTANT  CURVE  FIT 

Vincenti  and  Kruger  (Ref.  4,  p.  168  and  pp.  214-215)  recommend  an  Arrhenius-type  curve 
fit  for  the  equilibrium  constant  Kc,  specifically 

Kc  =  A  TB  exp(C/T)  (B-l) 

where  A,  B,  and  C  are  curve  fit  constants  to  be  determined  in  a  least-squares  sense  from 
the  equilibrium  constant  versus  temperature  data  for  a  given  chemical  reaction.  Table  B-l 
provides  the  A,  B,  and  C  curve  fit  constants  for  the  nine  chemical  reactions  of  present  interest 
based  upon  the  JANNAF  data  curve  fit  given  in  Ref.  13.  In  fitting  the  JANNAF  data  of 
Ref.  13,  the  fits  were  performed  over  a  low  (300-1 ,000  K)  and  high  (1 ,000-5,000  K)  temperature 
range  for  each  reaction  using  a  constrained  least-squares  curve  fit  which  forced  continuity 
of  the  fit  value  at  the  low-high  temperature  match  point,  namely  1,000  K.  All  calculations 
were  performed  using  double  precision  arithmetic  (16  digits)  with  all  JANNAF  data  input 
from  Ref.  13  also  performed  in  double  precision. 

Table  B-2  presents  the  standard  error  of  estimate  for  each  reaction’s  curve  fit  in  terms 
of  LN  (Kc)  where 


LN(KC)  =  LN(A)  +  B  LN(T)  +  C/T  (B-2) 

from  Eq.  (B-l).  Note  that  the  standard  error  of  estimate  is  of  the  order  0.01  to  0.00 1  for 
values  of  LN  (Kc)  in  the  range  -  10  to  —600,  which  clearly  shows  the  high  accuracy  of  the 
Arrhenius-type  curve  fit  relative  to  the  equilibrium  constant  data. 

Figures  B-l ,  B-2,  and  B-3  illustrate  graphically  the  behavior  of  the  equilibrium  constant 
data  for  each  reaction  when  plotted  in  terms  of  Eq.  (B-2)  variables,  namely  LN  (Kc)  versus 
1/T.  Figure  B-l  shows  that  the  H2  =  H  +  H  and  OH  =  O  +  H  reactions  have  nearly  the 
same  equilibrium  constant  variation  with  temperature;  similarly,  for  the  HzO  =  H  +  H 
+  O  and  NO2  =  N  4-  O  +  O  reactions  on  Fig.  B-2. 

For  high-temperature  air,  the  pioneering  work  by  Wray  (Ref.  22)  performed  in  1961 
provided  the  equilibrium  constant  curve  fits  for  the  reactions 


02-20 


N2  — 2N 


NO-N  +  O 
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which  are  used  in  the  present  Fig.  1  algorithm  pseudo-coding  for  high-temperature  ionized 
air.  Figure  B-4  shows  a  comparison  of  the  Wray  curve  fits  relative  to  the  present  JANNAF 
curve  fits  over  the  high-temperature  range  from  2,500  to  5,000  K.  Agreement  is  essentially 
exact  over  this  temperature  range. 
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02  =  0  +  0  H2  =  H  +  H  . OH  =  O  +  H 


Figure  B-l.  Chemical  equilibrium  constant  variation  with  temperature, 
Oj,  Hj,  and  OH. 
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N2  =  N  +  N  H20  =  H  +  H  +  O  . N02  =  N  +  O  +  O 


Figure  B-2.  Chemical  equilibrium  constant  variation  with  temperature, 
N2,  H20,  and  N02. 
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CO  =  C  +  O  - C02  = 


Figure  B-3.  Chemical  equilibriun 


i 


i 
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■e-  Wray  JANNAF 


1/T(1/°K) 

Figure  B-4.  Comparison  of  Wray  and  JANNAF  chemical  equilibrium  constant  curve  fit 
variation  with  temperature. 
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Table  B-l.  Equilibrium  Constant  Curve  Fit  Coefficients* 


■ 

Range 

A 

B 

C 

1 

L 

4.062702994580318 

0.2502842621822995 

-59606.12559562957 

H 

3047.358840886025 

-0.5810944837674182 

-60483.34670079666 

2 

L 

1.368670697255243 

0.3649632776333878 

-113312.4168571282 

H 

290.6847098177596 

-0.307468905987193 

-  114025.8191245152 

3 

L 

0.09447126526609221 

0.4675387830430657 

-52005.88558611053 

H 

21.6209968596366 

-0.1991760350848257 

-52833.50726295615 

4 

L 

0.03927314917445225 

0.4790269258033703 

-51064.14698292007 

H 

9.437740626012161 

-0.1982331074650644 

-51867.73124463861 

5 

L 

0.003362325850839724 

1.169921035885374 

-110586.4571538467 

H 

13843.73806527046 

-0.724239179515325 

-112732.7725377666 

6 

L 

2.802867695064791 

0.3794166065314485 

- 129073.3747337088 

H 

633.4027082089172 

-0.3045586535585327 

-129769.1043600149 

7 

L 

30136.53402335615 

-0.02686298150640948 

-  192979.4188364627 

H 

16648911484.32225 

-1.701127984990677 

-194636.1231254872 

8 

L 

0.6117641964723212 

0.2857408752663698 

-75604.48446764979 

H 

127.867751455553 

-0.3883093860925471 

-76290.71512510274 

9 

L 

217.6601569189995 

0.1482325072416944 

-112264.619822127 

H 

12708271.17655437 

-1.246715227962329 

- 1 13603.4909203498 

*  Kc  =  A(Tb)  exp(C/T) 

L:  Low-Temperature  Range  (300-1,000  K) 

H:  High  Temperature  Range  (1,000-5,000  K) 

r  Chemical  Reaction 

1  Oz  =  O  +  O 

2  Nj  =  N  +  N 

3  H2  =  H  +  H 

4  OH  =  O  +  H 

5  H20  =  H  +  H  +  O 

6  CO  =  C  +  O 

7  C02  =  C  +  O  +  O 

8  NO  =  N  +  O 

9  N02  =  N  +  O  +  O 


30 


AEDC-TR-91-1 4 


Table  B-2.  Equilibrium  Constant  Curve  Fit  Statistics* 


r 

Range 

Standard  Error  of 
Estimate  in  LN(KC) 

Maximum  LN(KC) 

Minimum  LN(KC) 

1 

L 

0.008446705 

-195.8516 

-56.47537 

H 

0.008945615 

-56.47537 

-9.038786 

2 

L 

0.004053238 

-375.3098 

-110.4775 

H 

0.004279492 

-110.4775 

- 19.75358 

3 

L 

0.0008926589 

-173.0451 

-51.1357 

H 

0.01417387 

-51.1357 

-9.214879 

wm 

L 

0.001472959 

- 170.7178 

-50.99236 

u 

H 

0.0103901 

-50.99236 

-9.833852 

5 

L 

0.009161721 

-367.6373 

-108.2001 

H 

0.02416324 

- 108.2001 

-19.21817 

6 

L 

0.005646218 

-427.0459 

-125.4218 

H 

0.002992297 

- 125.4218 

-22.09832 

n 

L 

0.02045509 

-633.0897 

-182.8515 

■ 

H 

0.007943118 

-182.8515 

-29.88921 

8 

L 

0.006034852 

-  250.8724 

-74.12206 

H 

0.003529303 

-74.12206 

-13.71652 

9 

L 

0.01920225 

-367.9731 

- 105.8577 

H 

0.004792126 

-105.8577 

-16.98331 

*  LN(Kc)  —  LN(A)  +  B  LN(T)  +  C/T  Definition  of  Standard 

L:  Low-Temperature  Range  (300-1,000  K)  Error  Qf  Estimate  (S): 

Maximum  at  300  K  and  Minimum  at  1 ,000  K 


H:  High-Temperature  Range  (1,000-5,000  K) 
Maximum  at  1,000  K  and  Minimum  at  5,000  K 


N 


X.fyp1  -  y t**? 


r  Chemical  Reaction 

1  02  =  O  +  O 

2  N2  =  N  +  N 

3  H2  =  H  +  H 

4  OH  =  O  +  H 

5  H2O  =  H  +  H  +  O 

6  CO  =  C  +  O 

7  C02  =  C  +  O  +  O 

8  NO  =  N  +  O 

9  N02  =  N  +  O  +  O 


s  = 


N  -  F 


N  =  Total  Number  of  Data  Points 
F  =  Total  Number  of  Functions  in  Curve  Fit 
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